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Abstract 

We describe and discuss a recently proposed quantum Monte Carlo algo- 
rithm to compute the ground-state properties of various systems of interacting 
fermions. In this method, the ground state is projected from an initial wave 
function by a branching random walk in an over-complete basis of Slater deter- 
minants. By constraining the determinants according to a trial wave function 
\iPt), we remove the exponential decay of signal-to-noise ratio characteristic 
of the sign problem. The method is variational and is exact if |V't) is exact. 
We illustrate the method by describing in detail its implementation for the 
two-dimensional one-band Hubbard model. We show results for lattice sizes 
up to 16 X 16 and for various electron fillings and interaction strengths. Be- 
sides highly accurate estimates of the ground-state energy, we find that the 
method also yields reliable estimates of other ground-state observables, such 
as superconducting pairing correlation functions. We conclude by discussing 
possible extensions of the algorithm. 

PACS numbers: 02.70.-c, Tl.lO.+x, 71.20.Ad, 71.45.Nt 
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I. INTRODUCTION 



We describe a new quantum Monte Carlo (QMC) algorithm to compute the ground-state 
(T = K) properties of systems of interacting fermions. Our method, which is approximate, 
removes the exponential scaling of computation time with system size that is characteristic 
of the infamous fermion "sign problem" in QMC simulations @-|5|. Here we discuss the 
general concepts of the algorithm and then describe details for its implementation using the 
Hubbard model as an example. The test results we present will show that the algorithm 
makes it possible to compute accurately, in computation times scaling algebraically with 
system size, general ground-state properties, such as superconducting pairing correlation 
functions. A brief description of the basic CPMC algorithm and some of the results on the 
Hubbard model were published earlier The algorithm, as it will be detailed here, can 
also be directly applied to study many other lattice models of electron correlations, such as 
the extended Hubbard model, the Anderson lattice model, etc., where computer simulations 
with existing QMC algorithms are often difficult and sometimes impossible. Application 
of the method to more general problem classes, such as atoms, molecules, and nuclei, is 
currently under study. 

The new algorithm, called the constrained path Monte Carlo (CPMC) method, has two 
main ingredients: The first is casting the projection of the ground state from an arbitrary 
initial state as importance-sampled branching random walks in a space of Slater determi- 
nants. The second ingredient, used only to deal with the sign problem, is constraining the 
paths of the random walks so that any Slater determinant generated maintains a positive 
overlap with a known trial wave function \iPt)- 

The first of the two ingredients is an exact procedure. As we will illustrate, it combines 
important advantages of two existing methods, the Green's function Monte Carlo (GFMC) 
I^J^Q] and the auxiliary-field quantum Monte Carlo (AFQMC) p|-p!2| methods. For exam- 
ple, our method shares with the latter the ease of computing expectation values of certain 
correlation functions, which are crucial to probe physical properties but which are often hard 
to compute accurately by the standard GFMC methods. On the other hand, it shares the 
GFMC concept of importance sampling with a trial wave function {ipr), which greatly im- 
proves its efficiency over the AFQMC method. In addition, the realization of the projection 
by open-ended random walks along the imaginary-time direction, in contrast to the standard 
AFQMC formulation, makes it practical and easy to implement the second ingredient, the 
constrained path approximation, and hence to eliminate the exponential scaling due to the 
sign problem. 

The constrained path approximation ensures that the Monte Carlo representation of 
the projected ground-state has no asymptotic signal-to-noise ratio decay in imaginary time. 
The resulting method is variational, with the computed ground-state energy being an upper 
bound, and becomes exact if \iPt) is exact. The constrained-path approximation builds 
upon the positive projection technique of Fahy and Hamann |13|, but can also be viewed 
as a generalization of the fixed-node pi|-p^ approximation in the GFMC method. Because 
of the different bases in which the approximations are applied, the effect of the constrained 
path approximation is expected to be different from that of the fixed-node approximation. 

In Section II, we will summarize the Green's function Monte Carlo and the auxiliary-field 
quantum Monte Carlo methods for ground-state calculations. Here, we will establish the 



2 



necessary concepts and formalisms from these existing approaches that are integral parts of 
our new method. In Section III, we describe the CPMC method in general terms, focusing 
on the concept of the importance-sampled random walks in Slater-determinant space, the 
nature and consequences of the constrained-path approximation, and the computation of 
expectation values. Implementations issues are discussed in Section IV in the context of 
the one-band Hubbard model. In Section V, we report results for this model that illustrate 
the accuracy and performance characteristics of the new method. Finally, in Section VI, we 
summarize and discuss several simple extensions of the CPMC method. 



II. BACKGROUND 

In this section, we summarize the AFQMC method and also sketch a particular GFMC 
method, namely the diffusion Monte Carlo (DMC) method |]14| , |T5| , that is most analogous 
to our new method. In discussing the DMC method, the approach we use is not stan- 
dard, but rather it is one designed to provide the necessary groundwork for the descrip- 
tion of CPMC. Both the AFQMC and the GFMC methods contain elements important to 
the CPMC method. For example, the basic techniques of the AFQMC method, such as 
Hubbard-Stratonovich transformation, imaginary-time propagation of Slater determinants, 
and matrix multiplication stabilization p,|lT],|T^, are shared by the CPMC method, on the 



other hand the random walk realization of the propagation []7|,[T8|, importance sampling in 
the random walks by use of a known trial function 0,0, and the fixed-node approximation 



IJ] are all GFMC concepts of much relevance. 

Most ground-state quantum Monte Carlo methods are based on 



|7/>o) oc jim e-^^|7/>r); (1) 



that is, the ground state lipo) can be projected from any known trial state \iPt) that satisfies 
(V'tIV'o) 7^ 0. In a numerical method, the limit can be obtained iteratively by 

|^(n+l)^ ^g-ArH|^(n)^^ (2) 

where \ip^^^) = lipr)- With a small At, the first-order Trotter approximation can be used: 

^-ArH ^ ^-ArK^-ArV_ (3) 

Typically, K and V are the kinetic and potential energy operators. More generally, they are 
the one- and two-body interaction operators. 



A. Auxiliary-field quantum Monte Carlo method 

In the AFQMC method, the operators and wave function are in a second quantized 
representation, defined in terms of fermion creation and destruction operators and c. The 
basis is one of Slater determinants: 

|0)^0t0t...0j,jO) (4) 
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where 



^\^Y.^]^,^. (5) 
j 

^ji are the elements of a matrix $ of dimension x A^^., where is the size of the basis 
and is the number of fermions with spin a. Each column of the matrix $ represents a 
single-particle orbital that is completely specified by a vector of dimension A^. One example 
of such a Slater determinant is the Hartree-Fock (HF) solution |(/)hf) = Ilcrl^HF)' where 
each |0hf) defined by a matrix $gp whose columns are the N„ lowest HF eigenstates. For 
any two real non-orthogonal Slater determinants, |0) and it can be shown that their 
overlap integral 



') = det($T$'), (6) 



and single-particle Green's function 

.ti 



" = - [^'i^^^r'<^%- (7) 

Now we consider the projection in this Slater-determinant basis. The trial wave 
function can be a linear combination of determinants, but without loss of generality, 
we assume it is a single determinant. A key point is that the projection of any Slater 
determinant by any operator of the form 



simply leads to another Slater determinant, i.e.. 



(8) 



ei:.,^lA^^.-.|0) = 0't0't...^,1^jO)^|0') (9) 

with <\)'\ = Y.j c] ^'ji and = e"*^*. 

The e~^'^^ part of (|^) has the single-particle form of (S). The e"^"^^ part, however, 
represents two-body interactions with V = 1/2 J^ijkiVijkic-CjCiCk- Following Hubbard, we 
rewrite V quadratic form: 



^ = E^"IE4^Sc.l =E^-pL (10) 



where the parameters and the matrix R" are defined by the elements Vijki and the 
number of a is at most A^^ but is often much smaller. With this quadratic form, a Hubbard- 
Stratonovich (HS) transformation of the two-body part of yields: 



Ar^^A^p^ =Y[r dxf-^ exp(x« v^-ArA„ ^ c\R%c,) , (11) 

„, J — CO \/ ZTT ■■ 



-oo V27r 



where is an auxiliary-field variable. Denoting the collection of such variables by x and 
defining B{x) = exp[— At clKijCj] Hq exp[xa\/—ATXa J2 ^iR'^jCj], we obtain: 
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e~^^^ = / dxP{x)B{x), (12) 



where P{x) = Y[a{^~^^°' / V'^^'^) is a probability density function and B{x) has the desired 
single-particle form of (P). The essence of the HS transformation is the conversion of an 
interacting system into many non-interacting ones living in fluctuating external auxiliary- 
fields, and the summation over all such auxiliary-field configurations recovers the correct 
many-body interactions. We note that different forms of this transformation exist and that 
they can affect the algorithm performance, possibly to a large degree. These issues are, 
however, not of concern here, as we will only be describing the general algorithm. 

With (0) and (0), the ground-state expectation (O) of some observable O can be com- 
puted by The denominator is 



2n 



^^(0)|g-nArHg-nArH|^(0)^ = / (t^^ | [fl t^x(')P(f« )S (f « )] (13) 

•' 1=1 

= f[l[dx^'^P{x'^''>)]det{^^]lB{x^'^)^T) (14) 

where B(x) is the N x N matrix associated with the single-particle operator B{x) and 
equations (|^), (H), and (0) have been applied. In the AFQMC method 0, is fixed and 
the many- dimensional integral in (0) is evaluated by a Monte Carlo (MC) method like the 
Metropolis algorithm. The MC process samples configurations {x^^\x^'^\ . . . ,5;*^^'^^} of the 
auxiliary-fields distributed according to the absolute value of the integrand. 

In the AFQMC method, the sign problem occurs because in general the determinant 
in (H^) is not always positive. In fact, its average sign approaches zero exponentially as n 
(or N) is increased The integral then becomes vanishingly small. Thus an exponential 
growth in computation time is required in its evaluation, since the MC samples, drawn 
from the absolute value of the integrand, become dominated by noise. This problem has 
remained largely uncontrolled, preventing general simulations at low temperatures or large 
system sizes. 

One attempt to control the sign problem was the positive projection approximation 



proposed by Fahy and Hamann [0. They used a known wave function \^pc) and imposed 
2n conditions 

(V'r|5(fW)5(f(^))---5(f('))|^/'e) > 0, / = l,2,...,n (15) 
(V^e|5(f('))B(f('+^))---5(f(2"))|V^T) > 0, / = 2n,2n- l,...,n (16) 

in sampling the auxiliary fields. The approximation is similar in spirit to that of the fixed- 
node approximation in the GFMC method. However, the constraint is non-local in imaginary 
time, as any change in x'-'^ affects the constraint conditions at a// times between / and n. Thus 
all auxiliary fields had to be updated simultaneously and only paths satisfying all constraint 
equations were accepted. The approach is hence computationally very intensive. In our 
CPMC method, we adopt the Fahy-Hamann concept of a constraining state but implement 
the constraint in the context of a random walk in the space of Slater determinants, which 
makes the procedure practical and straightforward. 
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B. Diffusion Monte Carlo method 



The DMC method |TP[ executes the iteration in as random walks in configuration 
space. When the fixed-point condition is reached, the random walks sample positions in 
configuration space from a distribution that represents the unknown amplitude of the ground- 
state wave function. 

We denote the configuration basis by \R), where R = {ri,r2, ■ ■ ■ ,rN^} is the electron 
coordinates in the continuous three-dimensional space. In this basis, the potential energy 
propagator e"^"^^ in @ is diagonal, but the kinetic energy propagator e~^'^^ , where K = 
— ■^J^i'^h is ^o^- order to write the latter in a more suitable form for a Monte Carlo 
treatment, we invoke a Hubbard- Stratonovich (HS) transformation: 



1,,2 
2-' 



Since e^'^^''^'|i?) displaces rj in \R) by V Arxi, the effect of e~^'^^ on any \R) can be viewed 
stochastically as "diffusing" it to \R + V Arx) , where each component Xi of the auxiliary 
field X is drawn from the normal distribution function P{xi) = e~2^'i /(27r)3/^. 

The wave function lip^"^^) can be expressed in terms of the amplitudes {R\ip^'^'') = il)^'^\R). 
In the random walk realization of the iteration, ip^"'\R) is represented by a finite ensemble 
of configurations {i?^"^}. At each stage, the Monte Carlo method provides the stochastic 
sampling of x^"^ and consequently the movement l-R^."'') l-R^"'' + V Atx^^'') = for 
each configuration in the ensemble. The factor e~^'^'^^^k ) translates into a weight (branching 
factor) for the configuration. As the iteration approaches the fixed-point condition, the 
weighted distribution of configurations represents ipo{R) = {R\ipo). 

The sign problem in the DMC method has a somewhat different character than the sign 
problem in the AFQMC method. The Pauli exclusion principle requires that the fermion 
wave function iPq{R) change sign if the positions of two electrons with the same spin are 
interchanged. Unlike the AFQMC method, the DMC method does not impose the anti- 
symmetric property in the projection process. Without additional mechanisms, the DMC 
method naturally produces points distributed according to the lowest eigenstate of the diffu- 
sion equation. This state is symmetric and bosonic-like. There has only been limited success 
in attempts to construct exact algorithms that yield asymptotically (in n) a non-vanishing, 
anti-symmetric Monte Carlo signal ||2yj2T|]. 



The fixed-node method |TJ-|T6[ is an approximate scheme to prevent the convergence 
to the bosonic-like ground state. Antisymmetry in ipo{R) implies that there are equivalent 
regions in configuration space that are separated by a nodal surface on which ipo{R) = 0. 
The exact nodal surface is in general unknown. In the fixed-node approximation, a trial 
nodal surface is assumed, based on a known trial wave function \iPt)- A solution which 
is everywhere positive is then sought in the region iPt{R) > by imposing the boundary 
condition that %jj^"\R) vanishes at ipriR) = 0. Unless ipriR) = happens to be the correct 
node, the resulting DMC solution for the ground state is approximate. The ground-state 



energy obtained is an upper bound |T^. 



An important feature of the DMC method is importance sampling. This technique is 
necessary to reduce the variance of the computed results to acceptable levels. For brevity 
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we will not discuss this technique here. Instead, we will postpone such a discussion until it 
is needed to complete our description of the CPMC algorithm. 



III. CONSTRAINED PATH MONTE CARLO METHOD 

We now describe the CPMC algorithm. It uses the Hubbard- Stratonovich-based for- 
malism of the AFQMC method, but a Monte Carlo sampling procedure similar to that of 
the DMC method. The iterative process (H) becomes an open-ended random walk in Slater 
determinant space. Within the framework of this random walk, we introduce importance 
sampling and the constrained path approximation. 

We remark that any antisymmetric wave function can be written as a linear combination 
of Slater determinants, i.e., 

I^)=EX40)I0), (18) 

where the sum is over each member of the Slater determinant basis. As introduced in Section 
II, we will always use to denote antisymmetric wave functions and |0) to denote a single 
Slater determinant. Contrary to the configuration space used in the DMC method, the 
Slater-determinant basis space of |0) is non- orthogonal and over- complete. 



A. Importance-sampled random walk formulation 

Using ([I2D , we write the iterative equation (Q) as 

= / rffP(f)5(x)|V^(")). (19) 



In the Monte Carlo realization of this iteration, we represent the wave function at each stage 
by a finite ensemble of Slater determinants, i.e., 

IV'^")) (X E (20) 

k 

Here k labels the Slater determinants and an overall normalization factor of the wave function 
has been omitted. The Slater determinants are referred to as random walkers as they are 
generated by the random walk. At any stage of the iteration, the sum will be over only 
part of the basis as the determinants in this sum are statistical samples whose distribution 
represents the linear coefficient x^^M iii The statistical accuracy of this representation 
increases algebraically as the number of independent samples is increased. In the remainder 
of the paper, equation (^) will serve as the definition of the Monte Carlo representation 
of a wave function in the CPMC method. We will start from an initial ensemble where 



One step of the iteration involves the propagation of each walker according to (plQl). Since 
the non-interacting operator B{x) operating on any Slater determinant leads to another 
Slater determinant, an analytical realization of this propagation for each walker would yield 
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a linear combination of many Slater determinants. In our random walk, this propagation is 
achieved stochastically by Monte Carlo (MC) sampling of x: 



t^'^)^ I dxP{x)B{x)\<pt\. (21) 



that is, for each random walker we choose an auxiliary-field configuration x from the proba- 
bility density function P[x) and propagate the walker to a new one via \(t)^k^^^) — ^{^)\'i^^k^)- 
We repeat this procedure for all walkers in the population. These operations accomplish 
one step of the random walk. The new population represents l?/'*^""'"^)) in the sense of (PUD. 
These steps are iterated indefinitely. After an equilibration phase, all walkers thereon are 
MC samples of the ground-state wave function and ground-state properties can be 
computed. 

In order to improve the efficiency of (|Typ and make it a practical algorithm, an importance 
sampling scheme is required. In the procedure just described, no information is contained 
in the sampling of x on the importance of the resulting determinant in representing |'?/'o), 
yet such information is clearly important. For example, the ground-state energy is given 
by Eq = {iPt\H\iPq) / {ipTl'ipo) ■ Hence, estimating £'0 requires estimating the denominator 
by J2(j,{'4'T\4') , in which |0) denotes random walkers after equilibration. Since these walkers 
are sampled with no knowledge of {iPt\4')j terms in the summation over can have large 
fluctuations that lead to large statistical errors in the MC estimate of the denominator, 
thereby in that of Eq. 

To introduce importance sampling, we iterate a modified equation with a modified wave 
function, without changing the underlying eigenvalue problem of (|19|). Specifically, for each 
Slater determinant we define an importance function 

0^(0) ^ (^^10), (22) 

which estimates its overlap with the ground-state wave function. We can then rewrite 
equation ([T9|) as 



= J dfp(f)E(f)|^/;(")), (23) 

where the modified wave function is 

|^(")) = ^Or(0)x^(.)(0)|0) (24) 

and the modified "probability density function" is 

^<^^' - i5!W^'">- '''' 

We note that P{x) is a function of both the future 10^"+-^-') and the current |0*^"')) positions 
in Slater-determinant space. It is trivially verified that equations (p!9| ) and p^ ) are identical. 

In the random walk, the ensemble of walkers { 10^"'') } now represents the modified wave 
function |'0*-"''), which is to say that their distribution represents the function x^mOt- The 
iterative relation for each walker is again given by (^1]), but with P{x) replaced by P{x). 
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The latter is in general not a normalized probability density function, and we denote the 
normalization constant for walker k by N{(j)^^^) and rewrite (21) as 



10^'^) - iV(4")) / dx^^BixM':^). (26) 

This iteration now forms the basis of the CPMC algorithm. It is convenient to associate 
a weight wl^^ with each walker, which can be initialized to unity. One step of the random 
walk is then as follows: For each walker 10^"''), («) sample a x from the probability density 
function P(f)/Ar(0i"^), (zz) propagate the walker by B{x) to generate a new walker, and 
(iii) compute a weight w^^'^^^ = tu^"''A^(0^"^) for the new walker. 

Steps (i) and (Hi) are sometimes difficult to implement. To ease their implementation, 
we apply the HS transformation of (^) and (|25|) to each component of x. This application is 
simple since both P{x) and B{x) can be decomposed into a product of independent factors 
corresponding to individual components Xa- Every step of the random walk then consists of 
successive sub-steps in which the Xa are sampled one by one, each according to {i) thru (iii). 
As we discuss in Section IV, such a decomposition is adequate to make the Hubbard-model 
application straightforward, since the HS transformation we use allows only two discrete 
values (±1) for each Xa, and thus the easy tabulation of P{xa)- For more general cases, 
however, it is often necessary to further approximate P{xa)- The following procedure can be 
adopted: Under the assumption of small At, the ratio of the overlap integrals is manipulated 
into the form of an exponential whose exponent is linear in Xa', P{xa) is then written as a 
shifted Gaussian times a normalization constant. The basic idea of this procedure is similar 
to that used in the DMC method. 

To better see the effect of importance sampling, we observe that if lipx) = I'^o); the 
normalization / P{x)dx is a constant. Therefore the weights of walkers remain a constant 
and the random walk has no fluctuation. Furthermore, we refer again to the estimator 
for Eq. With importance sampling, the denominator becomes the sum of weights w, while 
the numerator is Y.(f>{4'T\H\4>)w/ {ipT\(f>)y where again |0) denotes walkers after equilibration. 
As \iPt) approaches \ipo), all walkers contribute equally to the estimator and the variance 
approaches zero. We emphasize that different choices of importance functions only affect 
the efficiency of the calculation. 



B. Constrained path approximation 

Despite the advantages over the standard AFQMC method in terms of sampling effi- 
ciency, the random walk formulation still suffers from the sign problem. Here we will illus- 
trate the origin of the sign problem in this framework and then introduce the constrained 
path approximation to eliminate the exponential decay of the average sign. We will see that 
while such a fixed-node-like approximation has proved difficult to implement effectively in 



standard AFQMC method ||T3[, it is extremely simple to impose under our random walk 
formulation. 

The sign problem occurs because of the fundamental symmetry existing between the 



fermion ground-state and its negative —\ipo) For any ensemble of Slater deter- 



minants {\(f))} which gives a Monte Carlo representation of the ground-state wave function. 
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this symmetry implies that there exists another ensemble { — 10)} which is also a correct 
representation. In other words, the Slater determinant space can be divided into two degen- 
erate halves (+ and — ) whose bounding surface M is defined by (V^ol^) = and is in general 
unknown. 

In some special cases, such as the particle-hole symmetric, half-filled one-band Hubbard 
model, symmetry prohibits any crossing of M in the random walk. The calculation is then 
free of the sign problem ||2^ . In more general cases, walkers can cross M in their propagation 
by e~^'^^ . The sign problem then invariably occurs. Once a random walker reaches A/", it 
will make no further contribution to the representation of the ground state: 

{^(1)) = ^ = for any r. (27) 

Paths that result from such a walker have equal probability of being in either half of the Slater 
determinant space. Computed analytically, they would cancel, but without any knowledge 
of A/", they continue to be sampled in the random walk and become Monte Carlo noise. 
At sufficiently large n, the Monte Carlo representation of the ground-state wave function 
consists of an equal mixture of the -|- and — ensembles, regardless of where the random walks 
originated. The Monte Carlo signal is therefore lost. The decay of the signal-to- noise ratio, 
i.e. the decay of the average sign of {iPt\4>), occurs at an exponential rate with imaginary 
time. 

In this regard, the fermion sign problem appears very similar in either the DMC, 
AFQMC, or CPMC algorithms. The difference between the algorithms is that in the DMC 
algorithm minus signs appear when particles interchange positions in configuration space 
while in the CPMC and AFQMC algorithms the orbitals must interchange. The orbitals are 
an extended quantity and hence, at least for systems near a mean-field solution, the fermion 
sign problem is reduced. 

To eliminate the decay of the signal-to-noise ratio, we impose the constrained-path ap- 
proximation. It requires that each random walker at each step have a positive overlap with 
the trial wave function lipr)'- 

(^t|01"^) > 0. (28) 

This yields an approximate solution to the ground-state wave function, \iPq) = J2<j) 10); in 
which all Slater determinants |0) satisfy (|28|) . From (^Tf), it follows that this approximation 
becomes exact for an exact trial wave function \iPt) = li'o)- 

As a consequence of the constrained-path approximation, the ground-state energy Eq, 
computed by the estimator discussed in Section III A, is an upper bound to the true value 
Eq. To see this, we introduce an anti-symmetrization operator in the Slater-determinant 
space that extends any wave function defined in half the space by 10) into the whole 
space by ^^4, 10) ~ Z]-</, I ~ 0)- Since A^\iPq) is an eigenfunction of the modified Hamiltonian 
H'^ = H + V^, where is cxd at A/" and elsewhere, we have H'^^A^I'iPq)) = H (AtplipQ)) = 
E'q{A(I,\iPq)) . Both AfplipQ) and A^lipo) reside in the same Slater-determinant space and both 
are antisymmetric functions. Thus Eq > Eq. On the other hand, we recall that 
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(29) 



Therefore = E'q and E^>Eq. 

To implement the constrained-path approximation in the random walk, we redefine the 
importance function by (|22|): 

Ot(0) =max{(V^T|0),O}. (30) 
This prevents walkers from crossing the trial nodal surface M and entering the "— " half- 



space as defined by iV'r)- In the limit Ar — > 0, (|30D ensures that the walker distribution 
vanishes smoothly at M and the constrained-path approximation is properly imposed. With 
a finite Ar, however, P{x) has a discontinuity at M and the distribution does not vanish. 
We have found this effect to be very small for reasonably small imaginary-time steps Ar . 
Nonetheless, we correct for it by modifying P{x) near M so that it approaches zero smoothly 
at M. As we discuss in Section IV, the procedure is analogous to the mirror-correction ITBIJll 
used in the DMC method. 



C. Computing expectation values 

After the random walk has equilibrated, the distribution of random walkers represents 
the ground-state wave function \iIjq) under the constrained-path approximation. Various 
expectation values can then be computed from a population of these walkers and their 
weights. For example, the ground-state energy is 

_ {^t\H\xIj^) _ T.kMMH\<Pk)im(t^k) 

where terms in the numerator l(ipT\H\(f)k) / {'^T\'Pk) are given by combinations of elements of 
the Green's function as defined in (|^). 

An estimator similar to (pT|), namely {iI)t\0\iI)q) / (^/'tIV'o)) is easily obtained for any other 
operator O. In the GFMC method, this type of estimator is referred to as the mixed 
estimator. We recall that the true expectation value of O with respect to is 



For the energy, the mixed estimator is equivalent to the true estimator (|32D , but this is not 
true for any O that does not commute with H. In this case, it is sometimes possible to 
improve the mixed estimator by the following linear extrapolation 0: 

cxtrap ~ mixed (O)var, (33) 

where the variational estimate (O)var = ('^t|C'|V't)/(V^t|'^t)- 

Even a good trial wave function \iPt) with a good variational energy can sometimes fail 
to give a reasonable estimate for certain correlation functions. In such cases, ( |33| ) will not 



be effective. It is then imperative to compute (p2D. To do this, we devised a scheme called 
hack-propagation (BP), the essence of which comes from the forward walking (FW) technique 
1|] in the GFMC method: 
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r^oo {iPt exp{-rHc)\tp^) 

A subtle distinction, however, exists between back-propagation and forward walking. In 
back-propagation, ('?/^'rexp(— r/7c)| = (V'tI exp(— rife) is restricted to "constrained" paths, 
i.e., those paths that do not violate the constraint in the original forward direction 
exp(— Ariic)|V'o)- The standard ground-state DMC method has no sense of "time," because 
movements in configuration space are exactly reversible. The CPMC method, however, has 
a time direction: a set of determinants which does not violate the constraint at any time 
when going from right to left may indeed violate it any even number of times when evaluated 
from left to right. 

Because of this sense of direction, expression (0) may not yield (^2]). However, since 
\ipo) is itself approximate, this issue is not crucial. What is crucial is that (0)bp remains 
exact for an exact trial wave function. To demonstrate that it does, we will use perturbation 
theory and start by considering a Hamiltonian 

H' = H + AO, (35) 

where O is the operator whose expectation value we seek. We then apply the CPMC method 
to the new Hamiltonian H' with a new constraint governed by {iP't), where 

\ij'^) = \ijT) + X\SiJT), (36) 

and \6i/jt) is orthogonal to lipr), i-e., (V'tI'^t) = (V'tI^t), which is the standard boundary 
condition of perturbation theory. In the limit of small A, the constraint becomes identical to 
that of the original Hamiltonian H. To first order in A, we thus simply regain the previous 
expression (0) for (0)bp- The term proportional to \6iPt) does not contribute because it is 
orthogonal to the true ground state. 

Numerically we compared {0)bp with (O)'^ in several simple cases and observed rea- 
sonable agreement. In principle we can compute (O)'^ by creating a separate walk for the 
left-hand wave function, which propagates from with an appropriate sense of direc- 
tion, and then matching it with populations in the regular (right-hand) walk. In practice, 
however, it is difficult for such a scheme to yield accurate results, due to a lack of proper 



importance sampling |25 



Neither of the above measurement procedures is free from bias in the long imaginary- 
time limit. Since we are dealing with a branching random walk, there is necessarily a bias 
that arises from finite population sizes. However, this bias can be greatly reduced by taking 
a relatively large population size (a few hundred to a few thousand). The convergence and 
amount of bias will depend upon \iI)t) and the low-energy excitations of the system. We have 
found the back-propagation estimate of ground-state observables to be both statistically and 
physically accurate in our Hubbard model calculations. 

It is relatively simple to implement the back-propagation scheme on top of a regular 
CPMC calculation. We choose an iteration n and store the entire population { l^^"^) }. As 
the random walk proceeds, we keep track of the following two items for each new walker: 
(1) the sampled auxiliary- field variables that led to the new walker from its parent walker 
and (2) an integer that labels the parent. After an additional m iterations, we carry out the 
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back-propagation: For each walker / in the (n + m)**^ (current) population, we initiate a de- 
terminant {iPt\ and act on it with the corresponding propagators, but taken in reverse order. 
The m successive propagators are constructed from the stored items, with exp(— Ari^/2) in- 
serted where necessary. The resulting determinants | are combined with its parent from 
iteration n to compute ( (9) bp, in a way similar to the mixed estimator (0). The weights are 
given correctly by w\^^^\ due to importance sampling in the regular walk. Starting from 
another iteration n', this process can be repeated and the results accumulated. 



IV. IMPLEMENTATION ISSUES: THE HUBBARD MODEL 

The one-band Hubbard model is a simple paradigm of a system of interacting electrons. 
Its Hamiltonian is given by 

H = K + V = -tY^ {clc,, + c],c^^) + m^m^ , (37) 

{ij)cF i 

where t is the overlap integral, [/ > is the on-site Coulomb repulsion, rii^ = cl^Ci^, and ( ) 
indicates near-neighbors. We will take t = 1 and assume a two-dimensional square lattice 
of size N = L X L, with periodic boundary conditions. 

The physics of this model is rich, containing magnetism, a metal-insulator transition, and 
heavy Fermion behavior. Originally the model was proposed for ferromagnetism; today's 
interest focuses on the extent to which it might exhibit superconductivity away from the 
half-filled case of = = N/2. It is in the electron and hole-doped regions around half- 
filling that existing QMC methods experience a debilitating sign problem that restricts the 
simulations to small lattice sizes. In another paper, we will detail our study of the physical 
properties of this model . 



To illustrate the CPMC algorithm in more detail, we now describe our implementation 
of it for the Hubbard Hamiltonian. For this and related lattice models, we use the discrete 
version of the Hubbard-Stratonovich transformation 



J2 p(xi)e^"'("'T— a)^ (38) 



Xi = ±l 



where cosh(7) = exp(Art//2) and the probability density function p(xj) = 1/2 allows only 
Xi = ±1. Here we label components of the auxiliary field x by i, instead of by a, because of 
their one-to-one correspondence with lattice sites. 



A. Specific issues 

Each Slater determinant |0) = and similarly B{x) = B\x)B^{x). Any overlap 

integral between two Slater determinants involves the product of overlaps of individual spin 
determinants, e.g., {iPt\4>) = Yiai'^Tl^'') ■ Reflective of the interaction, the t and | spin de- 
terminants share auxiliary fields. Aside from this connection, they propagate independently. 
With these details about electron spin taken into consideration, all our previous discussions 
directly apply. 
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To reduce errors associated with the first-order Trotter approximation (0), we use the 
second order symmetric form, exp(— Arif) ^ exp(— Ari^/2) exp(— ArV^) exp(— Ari^/2). 
Thus, 

B{x) = Bk/2Bv{x)Bk/2, (39) 

where By{x) is the auxiliary- field-dependent propagator from the HS transformation. From 
i^), B^{x) = UibU^i), where 

with s(t) = 1 and = —1, and — e~^^^'^^'^. The probability density function is 

P{x)=UiP{x^). 

For a walker 10), we now describe one sub-step of the random walk in which the z*^ 
component Xi of the HS field is sampled and then bv{xi) is applied to the walker. If bv{xi)\(f)) 
is denoted by we have 

Thus the probability for picking one of the two possible values of Xi is p{xi)/[p{+l) 1)]. 
Once an Xi is chosen, |0) is propagated to obtain a new walker and the new weight is 
w [p{+l) 1)]. Since hy{xi) modifies only one row of the matrix it is straightforward 
and inexpensive to compute the resulting ratio of determinants for each of the two possible 
values of Xj, provided the inverse of the overlap matrix is known. As one moves 

to the next i, the overlap matrix can be efficiently updated by procedures that are almost 
identical to ones used in the AFQMC method |]9|,[TT||. 



As we mentioned in Section III B, with a finite Ar, equation (pOD causes p{xi) to be 
discontinuous at A/". To correct for this, we include a simple "mirror correction" procedure 
in the above sub-step: If walker |0) is close to Af, where we define "close" as when one 
application of 6y(xj) (with either Xi) would lead to a determinant with a negative 
overlap with lipr), we modify its weight to w/[l — {4't\(P') /Ot{(P)]- After the new walker 
has been accepted, we again check whether it is close to J\f. This is done by applying bv{xi) 
one additional time with the chosen Xj. Similar to w above, the weight w' is modified if the 
overlap {4'T\bv{xi)\(f)') is negative. We note that there is essentially no computational cost 
in computing this overlap. 

B. Algorithm outline 

The basic steps of the algorithm are: 

1. For each walker, specify its initial state to be some appropriate and assign its 
weight w and overlap Ot each a value of unity. 

2. If the weight of the walker is non-zero, propagate it via Bk/2'- 
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(a) Perform the matrix-matrix multiplications 

= 6^/2$ (42) 
and compute the new importance function 

O't = Ot(0')- (43) 

(b) If O'rp 7^ 0, update the walker, weight, and Or'- 

w^wO't/Ot, Ot^O't- (44) 

3. If the walker's weight is still non-zero, propagate it via Bv{x) = Jlj byi^j)'- 

(a) Compute the inverse of the overlap matrix 

Oi,. = (*?$)-!_ (45) 

(b) For each auxiliary field Xi {i = 1, N) 

i. compute p{xi) in which the ratio with Ot can be expressed in terms of Xi 
and the Green's function Gu as defined in (|^); 

ii. sample Xi and update the weight according to the discussion in Section IV 
A; 

iii. if the weight of the walker is still not zero, propagate the walker by by(xi) 
and then update Ot and Oinv; 

iv. apply the mirror correction in steps i or iii if necessary. 

4. Repeat step 2. 

5. Include an overall normalization factor in walker weight: w = w e^^^^, where Et is 
an estimate of Eq. 

6. Repeat 2 thru 5, which form one step of the random walk, for all walkers in the 
population. 

7. If the population of walkers has achieved a steady-state distribution, periodically make 
estimates of physical quantities. 

8. Periodically adjust the population of walkers. 

9. Periodically re-ort honor malize the columns of all the $ with non-zero weights. 

10. Cycle the process until an adequate number of measurements are collected. 

11. Compute final averages, estimate their statistical error, and stop. 



15 



We omitted the spin index a, but identical operations for both electron spin determinants 
are implied whenever a matrix manipulation is described. In presenting the above steps, 
we focused on illustrating the algorithm; the outline does not represent the most efficient 
implementation. 

We also ignored back-propagation. To implement it, an additional procedure is required 
in each random walk step in the measurement phase. As discussed at the end of Section III 
C, this procedure takes one of the following three possibilities: {i) storing all walkers (step 
6), (m) storing ancestry links (step 3) and auxiliary fields (step S.b.iii) for each walker, and 
copying them if necessary (step 8), or {in) back-propagating from (V'rl, computing (C)bp, 
and accumulating results (step 6). The length m of the back-propagation is pre-set as are 
the number of iterations in each measurement block, the frequency of measurements, etc. 

Additionally, we introduced several previously undiscussed steps and procedures. We 
now discuss these. 



C. Additional general algorithmic issues 

In this section, we discuss the use of £'t, the growth estimator, population control, and 
re-orthogonalization. Variants of the first three are present in all GFMC calculations, while 
the last one is adopted from the AFQMC method. We merely highlight these issues in the 
context of the new algorithm and point to references with more extensive discussions. While 
we will refer to the outline for the Hubbard implementation in Section IV B, the issues are 
general. 

As shown in step 5, a constant e^"^^^ multiplies the propagator e~^'^^ . If = Eq, the 
steady-state eigenvalue of equation (0) is unity. In other words, after equilibration the total 
weight of walkers remains a constant on the average. An inaccurate input of Et is only a 
minor concern, since the value of Et simply changes the total weight systematically by a 
constant factor. A statistical estimate of this factor is easily obtained from the random walk 
by observing the total weights through a number of iterations [^. On the average these 
totals scale as exp[— Ar(i?Q — Et)] and hence provide an estimate of Eq. This procedure for 
estimating the ground-state energy is often referred to as the growth estimator. 

In the random walk, one walker will eventually dominate all others and the random 
walk will spend most of its time sampling walkers which contribute little. To avoid such a 



loss of efficiency, a population control procedure |19| is needed (step 8). First, a branching 
(or birth/death) scheme is applied, in which walkers with large weights are replicated and 
ones with small weights are eliminated by some probability. There exist various ways to do 
this P, p!U| , pB| , with the guideline being that the process should not affect the distribution 
statistically. Branching allows the total number of walkers to fiuctuate and possibly become 
too large or too small. Thus as a second step, the population size is adjusted, if necessary, 
by rescaling the weights with an overall factor. Re-adjusting the population size introduces 
a bias and should only be done infrequently. If this is not possible due to a poor \iPt) or poor 
importance sampling, several calculations must be done with different (average) population 
sizes in order to extrapolate to the infinite population limit. 

The overall structure of CPMC resembles a typical GFMC calculation [Q. After equi- 
libration, we introduce an intermediate phase in which the growth estimator is computed. 
The length of this phase is a parameter and the outcome is used to adjust Et- The new 
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value of Et is then used in the next phase, which is divided into independent blocks. Be- 
cause of the correlated nature of successive steps, measurements are made only at suitable 
intervals. At an iteration n when measurements are taken, the average of a quantity O is 
= u^"') /d^'^\ In each block, we do not accumulate Instead we accumulate m*-"^ and 

ci'-"^ separately and compute a final average at the end of the block. The final result and 
statistical error estimate (step 11) are obtained from these bin averages. 

Step 9 is prompted by the fact that the repeated multiplication of B(x) leads to a 
numerical instability: relatively quickly numerical error grows to the point where 10^."^) 
represents an unfaithful propagation of \4>^k^)- This instability is well-known in the AFQMC 
method and is controlled |lTI| , p!7| by a numerical stabilization technique that requires the 
periodic re-orthonormalization of the single particle orbitals in 10^*''). In our calculation, we 
use the modified Gram-Schmidt procedure [0,0 which for each walker |0) factors the matrix 
$ as $ = QR, where Q is a matrix whose columns are a set of orthonormal vectors and R 
is a triangular matrix. After this factorization, $ is replaced by Q and the corresponding 
overlap Ot is replaced by Ot/ det(R). 



V. RESULTS 



We have presented a rather detailed discussion of the constrained path Monte Carlo 
method to make the CPMC algorithm as transparent as possible. In this section we present 
a variety of results from simulations of the Hubbard model that are chosen primarily for 
their importance in illustrating algorithmic issues. These results are for two-chain and 
two-dimensional lattices. First, we give results for the ground-state energy as a function 
of system size A^, filling fraction (n) = {N^ + Ni)/N, and trial wave function. Later, we 
examine results for other physical quantities such as pairing correlation functions. 



A. Ground-State Energy 



In Figure 1 we demonstrate the convergence and stability of the CPMC method. The 
energy (top) and overlap integral (bottom) are plotted as a function of imaginary time 
r = nAr. On the left, we demonstrate the initial convergence from the trial wave function 
to the ground state, and on the right, the asymptotic stability of the algorithm. This 
figure is constructed for an 8 x 8 lattice with (n) = 0.875, and U = 4. This parameterization 
generates a fairly difficult system for the AFQMC method because of the sign problem, {"^t) 
was an unrestricted Hartree-Fock wave function obtained with U = 0.4 [not 4) and yielded a 
variational energy of —53.05. The average number of walkers was 600. The behavior shown 
is quite typical of CPMC calculations. 

In the relaxation phase (r < 10), a trial energy of Et = —58.0 was used. This value 
is significantly higher than the true ground-state energy, so the overlap integral {ipT\ip^'^'') , 
and consequently the population size, grows quite rapidly. During this phase, population 
control is applied frequently (every 5 steps). The rapid fiuctuations in value of the overlap 
integral shown in the bottom left are due to re-adjustments of the total population size in 
each of these applications. For 10 < r < 20 (not shown), we used the same Et and applied 
population control with the same frequency, but computed the ground-state energy with the 
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growth estimator. This value was then used for the trial energy Et during the measurement 
phase (r > 20). In this phase, even though population control was applied every 10 steps 
and branching occurs, the total population size remained within the pre-set lower and upper 
bounds of 300 and 1200. Indeed, all vertical displacements in the overlap integral in the 
lower right of the figure correspond to the beginning of a block where we reset the population 
size to within 10% of the expected average of 600. 

Clearly, the values of both the overlap integral and the energy are completely stable as a 
function of n. In contrast, the standard AFQMC method has a bad sign problem. Indeed, 
for a 4 X 4 lattice at the same filling fraction and the same U, the average sign decays 
exponentially by roughly 5 orders of magnitude as the imaginary time increases from to 

20 i. 

The ground-state energy = —65.135 ± 0.008. (Extrapolation to At = will slightly 
reduce the average value.) The best AFQMC result available for this system is —65.02 ±0.06 
||27|| , which is slightly higher than our upper bound. The r-dependence of the energy after 
convergence is shown in the upper right of the figure. The inset demonstrates the short- 
term fluctuations in the mixed estimate of the energy; these fluctuations are also shown as 
the solid line in the main figure. Constructing local "binned" averages of the energy from 
blocks of length r = 10 yields the energy estimates shown as dotted lines in the inset. These 
binned averages are nearly statistically independent, and their root mean squared deviations 
divided by the square-root of the number of bins yielded the quoted statistical error. 

The computational requirements of the CPMC method are fairly modest when compared, 
for example, with the AFQMC method. Efficiency, however, can be dramatically affected 
by implementation issues. Two important issues are the accuracy of the trial wave function 
and of the imaginary-time propagator. To obtain the maximum efficiency, it is essential to 
use a propagator accurate to second order in the break-up between kinetic and potential 
terms. For example, figure 2 shows the energies obtained for first-order and second-order 
Trotter approximations for the propagators for a 4 x 4 lattice with N-^ = Ni = 6 spins and 
U = 4. The first-order propagator introduces significant systematic effects even for quite 
small time steps, however, the second order propagator permits significantly larger values 
of At to be used. We note that we also included a mirror correction to the propagation to 
make certain that there are no corrections of order less than (Ar)^. In this case, however, 
the differences between using second-order propagators with and without this correction are 
comparable to the error bars in the figure. 

Of course, while stability and efficiency are necessary for a useful simulation, accuracy 
is the principal concern. In Tables I and II, we compare our results for the Hubbard model 
to exact results for small systems and to other methods for large systems. Since the CPMC 
method depends upon approximate knowledge of the ground state of the system, we also 
included results for different trial wave functions. To date, we have used only free-electron 
and unrestricted Hartree-Fock (uHF) trial wave functions. 

In Table I, we see that the accuracy of the CPMC ground-state energy is always better 
than 5%, often much better, even when the trial wave function is very poor. The worst 
case is 3 X 3 with 4 | and 4 | spins, which is an open-shell case that corresponds to a 
very difficult filling fraction for the AFQMC method, and has a large f/ of 8, which makes 
single-determinant trial wave functions rather poor approximations. The Hartree-Fock wave 
function used as \iPt) had a very poor energy, —0.0025, compared to the exact energy of 
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—0.809. The CPMC method however was still able to obtain an energy of —0.766(2). Clearly 
a more accurate approximation to the ground state wave function would yield an even better 
result. 

In Table II we compare our results with available data from other numerical approaches, 
including stochastic diagonalization (SD) [Q, AFQMC, and density-matrix renormalization 
group (DMRG) |3ll] methods. The SD method uses Monte Carlo methods to attempt the 
construction of an efficient basis for approximating the ground state wave function of the 
system. Since an explicit basis is used, no sign problem occurs; however, an exponential 
growth in computing time occurs reflecting the increased effort in selecting members of the 
basis as system size increases. In contrast, the AFQMC method is in principle exact, but 
suffers from exponential growth in computing time as system size increases because of the 
sign problem. Finally, the DMRG method is a variational method that is very effective for 
one-dimensional and quasi-one-dimensional models. 

In Table II we focus filling fractions (n) around or greater than 80% as they are more 
interesting and also more difficult for QMC calculations. For the closed-shell 4x4 system, 
the results from the CPMC, SD, and AFQMC methods agree well. As we increase the lattice 
size, the SD results are comparatively poorer, presumably because of an insufficient number 
of states. The results from the AFQMC and CPMC methods continue to agree well up to 
fairly large lattice sizes. The worst case is an 8 x 8 lattice with 25 | and 25 J, spins, where 
the CPMC result lies approximately 0.4% ± 0.2% above the AFQMC result. For still larger 
lattice sizes, the sign problem limits the use of the AFQMC method to only closed-shell 
systems, for which the sign problem is much reduced. Even in these cases, we see that the 
error estimates in AFQMC are much larger than the corresponding statistical errors from 
the CPMC method. In the 12 x 12 case, the difference is roughly a factor of 30, i.e., about 
a factor of 900 more in CPU time. (Our calculations typically took tens of hours on an 
IBM RS6000 590.) Furthermore, the CPMC result is actually lower in energy in this case 
than the AFQMC results, but size of the error bars in AFQMC are similar to the difference 
between the two results. The CPMC result on a 16 x 16 lattice, a size far beyond the reach 
of the AFQMC method due to the sign problem, was obtained from a simulation comparable 
to that for the 12 x 12 system in terms of the numbers of iterations and walkers. For the 
two-chain system, we obtained excellent agreement with the DMRG results of Noack et al. 
Here the energy agreed to within less than 0.1%. 

B. Correlation Functions 

To effectively study ground-state properties, we need to accurately calculate pairing cor- 
relation functions, momentum distributions, and other ground-state expectation values. In 
previous sections, we discussed various ways of estimating a ground-state expectation value, 
in particular the back-propagation scheme. Here we will mostly benchmark results and 
further discuss related algorithmic behavior. We remark that for any simulation method 
correlation functions are in general much more difficult to compute than the energy. Thus 
there is limited data available from other methods with which we can benchmark our cor- 
relation functions. This means self-consistency checks (e.g., comparison of results obtained 
with different choices of IV't)) are crucial. 

In Table III, we show results on the simple closed-shell 4x4 system of 5 t 5 J, electrons 
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at f/ = 4 and with At = 0.05. The one-body density matrix is the expectation value of 
the Green's function elements: = (cgCi), where 1 = {lx,ly)- The spin density structure 
factor is 



where Si = ni| — riii is the spin at site 1. The charge density structure factor Sd is similar 
to with spin replaced by density, i.e., with the — sign in Si replaced by a + sign. The 
electron pairing correlation is defined as 



where a indicates the nature of pairing. The on-site s-wave pairing function has Ais(l) = 
cijCij, while in this case for ci-wave we used A2d(l) = ci-f J^s f{^)c\+si, where S is (±1, 0) and 
(0, ±1). For 6 along the a;-axis, f{6) is 1; otherwise it is —1. We average over different 
sites to improve statistics. 

In Table III we compare CPMC results from the back-propagation estimate with exact 
results. (The length of back-propagation was r = 6.) We see that the CPMC result is essen- 
tially exact for this system. The variational results suggest that the free-electron is not 
a very good trial wave function. Nonetheless the constrained-path error seems to be negligi- 
bly small. In fact, very limited branching occurs in the calculations, which indicates effective 
importance sampling with \iPt)- Furthermore, as U is increased to 8, the variational results 
become worse, but the CPMC results obtained with the same {iPt) remain accurate ||^. Also 
in this table are estimates by the mixed and extrapolation schemes discussed in III.C. While 
these schemes often improve the variational results considerably, their systematic errors are 
significant, particularly when compared to the high level of statistical accuracy that can be 
achieved with the CPMC method. 

In Table IV, we show expectation values for a system of 7 t 7 | electrons. This open- 
shell case has the worst sign problem for a 4 x 4 system. We show results from CPMC 
simulations with two different trial wave functions. Both are unrestricted Hartree-Fock wave 
functions, but {ipTi) was obtained with a U of 0.1, while \1pT2) with U = 4. The calculation 
with iV'ri) has much less fluctuation, even though \tpT2) has a lower variational energy. In 
fact, we found this trend to be rather general: free-electron-like wave functions tend to be 
better importance functions than unrestricted Hartree-Fock wave functions. We see that the 
two trial wave functions yield very different variational estimates, but their CPMC results 
are consistent and in reasonable agreement with exact results. For example, in the free- 
electron-like {ipTi), the momentum distribution is a step function, and the k = (1,0) state 
is completely occupied (n^ = 1), but with this trial wave function the CPMC method still 
gives the correct occupation of 0.92(1). 

In Figure 3 we show the (i-wave pairing correlation function D2dO) and static magnetic 
structure factor S'(k) obtained with two different trial wave functions for a 6 x 6 system. 
The expectation values obtained directly from these two different trial wave functions are 
shown in thick lines, and the corresponding back-propagation estimates {0)bp are shown 
with thin lines. While the two CPMC estimates do not agree exactly, they do demonstrate 
a much closer correspondence with each other than those obtained with the original wave 
functions. The case shown is comparatively easy because of the small size of the lattice and 
the relatively low value of electron filling. However, the overall trend is rather general. 



S{K,ky) = S{k) = l/iV^exp(zk-l)(soSi) 



(46) 



D,(/,,g = D(l) = (At(l)A,(0)) 



(47) 
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VI. SUMMARY AND DISCUSSION 



We described in detail the background, formalism, and implementation of the new con- 
strained path Monte Carlo algorithm. The CPMC method is a general quantum Monte 
Carlo algorithm for computing fermion ground-state properties. It introduces several new 
concepts, including importance-sampled random walks in a Slater-determinant space and 
the constrained path approximation within this framework. The algorithm combines advan- 
tages of the existing Green's function Monte Carlo and auxiliary-field quantum Monte Carlo 
methods, is free of any signal-to-noise ratio decay, and scales algebraically with system size. 
Together with data in Ref. we demonstrated that the method produces very accurate 
results for the Hubbard model, even with very simple choices of the trial wave function \'4)t)- 

Compared to the GFMC method, the current algorithm allows the random walk to take 
place in a basis other than that of configurations or occupation numbers. In this sense, the 
CPMC algorithm is a generalization of the GFMC algorithm. The CPMC method expresses 
the ground-state wave function (stochastically) in the form of (pOD, rather than ipo{R) = 
J2k ^(-R ~ ^k) as in the GFMC method. This is advantageous as it makes feasible the use of 
various techniques developed for one-electron calculations for atoms and solids. In addition, 
it makes our back-propagation scheme efficient and effective. Thus expectation values can 
be computed via (|^), while in the GFMC method the analogous forward walking technique 
has often been difficult and computations of some correlation functions have almost been 
impossible. If applications of the CPMC method to continuum systems are successful, the 
ability to compute expectation values such as forces will be very valuable. Such applications 
are under study. 

There is an obvious resemblance between the constrained path (CP) approximation in the 
CPMC method and the fixed-node (FN) approximation in the GFMC method. Both result 
in solutions to the Schrodinger equation that are consistent with some artificial boundary 
conditions. An important difference, however, is also evident. The FN approximation is in 
configuration space and requires the solution to have a pre-defined node: ip™{R) = where 
iPt{R) = 0. The CPMC method is in a Slater determinant space. The CP approximation on 
each individual Slater determinant translated into configuration space \R), is non-local 
and requires / ipT{R)<P{R)dR > 0. Thus the node, as well as the amplitude 0(-R), is allowed 
to vary. The systematic error in the CPMC method arises because the solution it yields, in 
the form of (0), has the artificial constraint XV'o(<^) > 0- 

In the CPMC method, each Slater determinant |0) analytically defines a continuous 
function (j){R) in contrast to walkers in the GFMC method which are delta-functions. It is 
thus easier to impose symmetries. One example is the trivial case of a non-interacting system: 
the CPMC method naturally yields the correct result, while standard GFMC still requires 
knowledge of the node. Another is the one-band Hubbard model at half-filling, where the 
CPMC method retains the exact nature of the AFQMC method, while the GFMC method 



retains the sign problem . 

Recently, ten Haaf and van Leeuwen presented data on the Hubbard model from 
standard GFMC simulations with the fixed-node approximation, which they and collabora- 
tors had earlier generalized to treat lattice fermion systems. For the 4x4 system in 
Table II (5 t 5 j f/ = 4), the CPMC result for the energy per site is E/N = -1.2239(3) 
(exact value [^: —1.2238). With an identical trial wave function, the fixed-node calcula- 
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tions of ten Haaf and van Leeuwen yielded —1.2186(4). Incorporating a Gutzwiller factor 
only slightly improved their FN result to —1.2201(4). Unfortunately, the rest of their re- 
sults are all FN energies computed for half-filled systems, for which the CPMC energies 
would be exact. Further comparisons away from half-filling would allow a more systematic 
understanding of the relative strengths of the CPMC method. 

It is worth noting that the CPMC algorithm provides a stochastic method closely linked 
with more traditional quantum chemistry approaches such the configuration interaction (CI) 
method. Similar to the CI method, the CPMC method produces a collection of determinants 
whose sum represents the ground-state wave function. The determinants, however, do not 
have to be orthogonal to each other. Furthermore, they are generated efficiently and sys- 
tematically by a Monte Carlo process that is guided by importance sampling. The drawback 
of the CPMC method is of course its variational nature due to the CP approximation. 

We are currently investigating several schemes for further improving the algorithm. These 
include a method analogous to released-node technique |TH| in the GFMC method. For the 



energy, this seems straightforward. For other expectation values, it involves evaluating 
(C')bp(ti, T2) defined as: 

(0)Bp(ri,r2) = (V-T exp[-r2if^exp[-riiJ]|0|exp[-riif]|^o). (48) 

In this expression H is the original Hamiltonian without the constraint, and indicates 
the Hamiltonian in the presence of the constraint. Hence, for a period of twice ri, we evolve 
the system without constraint and for another T2 we include the constraint. In the limit 
of zero Ti we obtain the approximation used to date, while finite ti improves the estimate, 
i.e., makes it exact, at the cost of increasing statistical error. The bookkeeping in such a 
calculation could be arranged to calculate directly the difference between the current (C)bp 
and the transient estimation, and hence to provide a stringent test on the accuracy of a 
given calculation. 

Other possibilities for improving estimates of expectation values include optimization 
techniques for improving the trial wave function. As mentioned earlier, the algorithm as 
described can be used with a multi-determinant |V't), with the computational cost increasing 
only linearly with the number of determinants. Thus it is desirable to have good trial 
wave functions in the form of linear combinations of Slater determinants. In addition, wave 
functions in this form that can be tuned systematically to yield different properties would be 
highly useful, since self-consistency checks with the CPMC method can then be carried out 
simply by changing parameters in \'4)t)- Yet other algorithmic topics include the development 
of interacting- walker PD[ and mirror potential P7| analogs. 
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TABLES 



TABLE L Hubbard model ground-state energies per site from CPMC simulations compared 
with exact results. The first column under "system" is the lattice size; the second, the numbers 
of electrons with | and J, spins. U is the on-site Coulomb repulsion. The trial wave function \iPt) 
used in CPMC is either a free-electron (free) or an unrestricted Hartree-Fock (uHF) wave function. 
E'var indicates the corresponding variational energy from this wave function. Statistical errors are 
in the last digit and are shown in parentheses. Exact results for the 4x4 systems are taken from 



[28-3C1. 



system 


u 










2x2 


2T U 


4 


uHF 


-1.5327 


-1.6038(6) 


-1.6046 


2x3 


2T 2i 


4 


free 


-1.267 


-1.3828(9) 


-1.4009 


2x3 


2T 2i 


8 


free 


-0.889 


-1.221(2) 


-1.244 


2x4 


2T 2i 


4 


uHF 


-1.333 


-1.3678(5) 


-1.374 


2x4 


3T 3i 


4 


uHF 


-1.438 


-1.5693(5) 


-1.569 


1x8 


3T 3i 


4 


free 


-0.645 


-0.8329(7) 


-0.834 


3x3 


4T 4i 


8 


uHF 


-0.0025 


-0.766(2) 


-0.809 


4x4 


2T 2i 


4 


uHF 


-2.8225 


-2.8813(3) 


-2.8825 


4x4 


4T 4i 


4 


uHF 


-1.025 


-1.095(1) 


-1.096 


4x4 


5T 5i 


8 


free 


-0.7188 


-1.0925(7) 


-1.0944 


4x4 


6T 6i 


4 


uHF 


-1.3117 


-1.4763(5) 


-1.478 


4x4 


7T 7i 


4 


uHF 


-0.8669 


-0.9831(6) 


-0.9838 


4x4 


7T 7i 


12 


uHF 


-0.474 


-0.606(5) 


-0.628 



TABLE IL Hubbard model ground-state energies from CPMC simulations compared with avail- 
able results from other approaches. The first two columns follow the same convention as the cor- 
responding ones in Table I. The interaction strength U is 4. The stochastic diagonalization (SD) 
results are from |3C|; the density-matrix renormalization group (DMRG) results on two-chains are 



from [32|. The statistical errors are in the last one or two digits, as indicated. 



system 








-E'afqmc 


4x4 5T 5i 


free 


-19.582(5) 


-19.58 


-19.58(1) 


6x6 13T 13i 


free 


-42.34(2) 


-40.77 


-42.32(7) 


6x6 141 14i 


uHF 


-40.17(2) 




-40.44(22) 


8x8 25T 251 


free 


-72.48(2) 


-67.00 


-72.80(6) 


8x8 27T 27i 


uHF 


-67.46(4) 




-67.55(19) 


10x10 411411 


free 


-109.55(3) 




-109.7(6) 


12x12 611611 


free 


-153.43(5) 




-151.4(1.4) 


16x16 lOlt lOli 


free 


-286.55(8) 






system 


W) 


-EcPMC 


-E'dmrg 




2x8 7T 7i 


free 


-13.067(4) 


-13.0664(2) 




2x16 141 141 


free 


-26.87(2) 


-26.867(3) 
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TABLE III. Comparisons of the computed expectation values and correlation functions by 
different estimators for a 4 x 4 lattice with 5 t 5 J, electrons and [/ = 4. The trial wave function 
is the free-electron wave function. The corresponding variational values from it are shown in the 
first row. The next row contains mixed estimates from CPMC simulations, while "extrp" shows 
the values extrapolated from the first two rows via (|^). The row labeled BP gives the CPMC 
result via the back-propagation scheme. In the last row, Dig is from AFQMC simulations |3C|, 
while the others are exact diagonalization results, with the first four from and the last from 
p3| . Ek indicates the kinetic energy, p{lx, ly) the one-body density matrix, S and Sd the spin and 
charge density structure factors, and Dig and the s- (on-site) and d-wave pairing correlations 
respectively. The Monte Carlo errors are shown in parentheses. 





Ek 


P(2,l) 


5(7r,7r) 


S'd(7r,7r) 




D2d{2,l) 


variational 


-24.0 


-0.0625 


0.625 


0.625 


0.003906 


0.03125 


mixed 


-24.0(0) 


-0.0625(0) 


0.6938(4) 


0.5572(1) 


0.000684(3) 


0.03095(2) 


extrap 


-24.0(0) 


-0.0625(0) 


0.763(1) 


0.4894(2) 


-0.002538(6) 


0.03065(4) 


BP 


-22.55(2) 


-0.0563(3) 


0.729(1) 


0.508(1) 


-0.000615(9) 


0.0246(2) 


exact 


-22.52 


-0.0560 


0.73 


0.506 


-0.00058(5) 


0.02453 



TABLE IV. Computed expectation values and correlation functions from CPMC for a 4 x 4 
lattice with 7 | 7 | electrons and C/ = 4, compared with exact results. Results are shown for 
two different trial wave functions \ipTi) and \'4'T2)- Exact diagonalization results are from |28|; 



numbers in parentheses indicate either the range of values due to the ground-state degeneracy or 
uncertainties in extracting numbers from a graph. Statistical errors on the last digit of the CPMC 
results are in parentheses. Symbols are the same as in Table III. rik is momentum distribution. 







Ek 


/'(1,0) 


p(2,2) 


^(Tr, 7r) 


5'd(7r,7r) 


nfc(7r/2,0) 


variational 


I^Tl) 


-24.0 


0.1875 


-0.0625 


1.654 


0.625 


1.0 




|'0T2) 


-21.88 


0.1706 


-0.0602 


4.39 


0.516 


0.941 


CPMC 


I^Tl) 


-21.44(2) 


0.168(1) 


-0.051(1) 


2.90(1) 


0.432(1) 


0.92(1) 




\'4^T2) 


-21.39(8) 


0.168(1) 


-0.049(1) 


2.92(2) 


0.430(1) 


0.92(1) 


exact 




-21.39(1) 


0.168(1) 


-0.051 


2.16(2) 


0.425 


0.93(1) 
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FIGURES 
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FIG. 1. Stability and convergence of the CPMC method. In the lower half of the figure is the 
overlap integral (V'tIV'^"^) ^ a function of the imaginary time r = nAr, where n labels the random 
walk step. In the presence of the sign problem, this overlap integral would decay exponentially with 
n. In the upper half of the figure is the corresponding estimate of the total energy as a function 
of n. The inset is an enlarged version of the portion between r = 37.5 and r = 117.5 in which the 
dotted line indicates the computed energy value from blocks of length r = 10. This calculation 
yielded a ground-state energy of Eq = —65.135 it 0.008, compared to an AFQMC result |27| of 
-65.02 ±0.06. 
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FIG. 2. Illustration of the Trotter approximation error. The computed energy per site is shown 
as a function of the time-step At for a 4 x 4 lattice with 6 t 6 J, electrons and U = 4. The 
first-order Trotter approximation leads to significant finite time-step error. With the second-order 
approximation, convergence to the At = limit is much more rapid. The right-triangle indicates 
the exact energy for this system. The variational nature of the CPMC energies is visible. The 
Monte Carlo error bars are indicated. Curves are to aid the eye. 
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FIG. 3. Sensitivity of CPMC results to the choice of trial wave function for some corre- 
lation functions. The upper figure is the d-wave electron pairing correlation D2d(l) and the lower 
curve is the magnetic structure factor <S'(k). The frcc-clectron wave function \ipTi) and the unre- 
stricted Hartree-Fock wave function \1pT2) (with U = 4) are used. The corresponding mean-field 
results for these correlation functions are also shown (thick lines). The computed ground-state 
energies from CPMC are -42.345(3) and -42.295(16) (cf table II). 



30 



